/* OS_Multi_Array.h
 *
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
#include <src/miscellaneous/OS_basic.h>
#include <boost/operators.hpp>
#include <iostream>
/*! \file OS_Multi_Array.h
    \brief Multidimensional array specialized for scientific computing.
	
	This file define classes acting basically as std::vector<std::vector<....<T>..>>
	where T is of integral type (double, int etc...).
	These classes are designed for scientific computing :
	- the memory is organized as a contiguous block of data. Thus the copy is very fast.
	- they are designed such that the operator [] can be overloaded. This property is important to handle boundary conditions.
	- A full set of algebra operations is provided and optimized for performances. It uses mainly iterator incrementations.
*/
#pragma warning(disable:4996)
#pragma warning (disable: 4715)
#pragma warning (disable: 4172)
 
#if !defined(_OS_multi_array_)
#define _OS_multi_array_

/*! 
	Open issues : 
Performance issues on reference / dereference.

class OS_data_array : 
declare
		T*						OS_data_array::first_;
		indice_list_type		OS_data_array::size_;
		capacity_list_type		OS_data_array::capacity_;
as private members of OS_data_array. This links to template friends definition problems.

class OS_data_array_ref : 
	declare 
	mutable			smart_ptr<data_type>	OS_data_array_ref::ptr_;
	mutable			T*						OS_data_array_ref::first_;
as private members of OS_data_array_ref.
*/
//////////////////////////////////////////////////////////////////////////////
/*! \brief optimized minus operator A-B
	perform the algebra operation A-B, where A and B are two multidimensional vectors.
*/
template<class A,class B>
inline A& OS_Array_minus(A& left, const B& right)
{
	typedef typename A::size_type size_type;
	typedef typename A::base_value_type base_value_type;
	for (int i = 1;i<= A::Dim;++i)
		OS_DYNAMIC_CHECK(left.size(i) == right.size(i),"incompatible size in OS_array template inner substraction"); 
	base_value_type * begin = left.first_;
	base_value_type * from_begin = right.first_;
	const base_value_type * end = left.first_+left.capacity();
	for (;begin != end; ++begin, ++from_begin)
		*begin -= *from_begin;
	return left;
};
/*! \brief optimized minus operator A-b
	perform the algebra operation A-b, where A is a multidimensional vectors and b is a scalar.
*/
template<class A>
inline A& OS_Array_minus(A& left, typename A::base_value_type from)
{
	typedef typename A::base_value_type base_value_type;
	base_value_type * begin = left.first_;
	base_value_type * end = left.first_+left.capacity();
	for (;begin != end; ++begin)
		*begin -= from;
	return left;
};

/*! \brief optimized minus operator A+B
	perform the algebra operation A+B, where A and B are two multidimensional vectors.
*/
template<class A,class B>
inline A& OS_Array_plus(A& left, const B& right)
{
	typedef typename A::size_type size_type;
	typedef typename A::base_value_type base_value_type;
	for (int i = 1;i<= A::Dim;++i)
		OS_DYNAMIC_CHECK(left.size(i) == right.size(i),"incompatible size in OS_array template inner addition"); 
	base_value_type * begin = left.first_;
	base_value_type * from_begin = right.first_;
	const base_value_type * end = left.first_+left.capacity();
	for (;begin != end; ++begin, ++from_begin)
		*begin += *from_begin;
	return left;
};
/*! \brief optimized minus operator A+b
	perform the algebra operation A+b, where A is a multidimensional vector and b is a scalar.
*/
template<class A>
inline A& OS_Array_plus(A& left, typename A::base_value_type from)
{
	typedef typename A::base_value_type base_value_type;
	base_value_type * begin = left.first_;
	base_value_type * end = left.first_+left.capacity();
	for (;begin != end; ++begin)
		*begin += from;
	return left;
};
/*! \brief external product A^B
	Perform the algebra operation A^B, where A and B are two multidimensional vectors.
	If \f$ A = \{a_i}_i\f$ and  \f$ B = \{b_i}_i\f$, then it corresponds to the outer product \f$ A^B = \{a_i \cdot b_i\}_i \f$.
*/
template<class A,class B>
inline A& OS_Array_xor(A& left, const B& right) {
	typedef typename A::size_type size_type;
	typedef typename A::base_value_type base_value_type;
	for (int i = 1;i<= A::Dim;++i)
		OS_DYNAMIC_CHECK(left.size(i) == right.size(i),"incompatible size in OS_array template inner addition"); 
	base_value_type * begin = left.first_;
	base_value_type * from_begin = right.first_;
	const base_value_type * end = left.first_+left.capacity();
	for (;begin != end; ++begin, ++from_begin)
		*begin *= *from_begin;
	return left;
};
/*! \brief external product A/B
	Perform the algebra operation A/B, where A and B are two multidimensional vectors.
	If \f$ A = \{a_i}_i\f$ and  \f$ B = \{b_i}_i\f$, then it corresponds to \f$ A/B = \{a_i / b_i\}_i \f$.
*/
template<class A,class B>
inline A& OS_Array_xor_not(A& left, const B& right) {
	typedef typename A::size_type size_type;
	typedef typename A::base_value_type base_value_type;
	for (int i = 1;i<= A::Dim;++i)
		OS_DYNAMIC_CHECK(left.size(i) == right.size(i),"incompatible size in OS_array template inner addition"); 
	base_value_type * begin = left.first_;
	base_value_type * from_begin = right.first_;
	const base_value_type * end = left.first_+left.capacity();
	for (;begin != end; ++begin, ++from_begin)
		*begin /= *from_begin;
	return left;
};
/*! \brief tensor product A^B
	Perform the algebra operation A^B, where A is a matrix and B is a vector.
	If \f$ A = \{a_{i,j}}_i\f$ and  \f$ B = \{b_i}_i\f$, then it returns the matrix \f$ \{ a_{i,j} \cdot b_i\}_{i,j} \f$.
*/
template<class A,class B>
inline A OS_Array_Tensor(const A& left, const B& right) {
	OS_STATIC_CHECK(A::Dim == 2*B::Dim);
	typedef typename A::const_iterator	A_const_iterator_type;
	typedef typename A::iterator		A_iterator_type;

	A result(left);
	A_iterator_type			A_iterator				= result.begin();
	A_const_iterator_type	A_end					= result.end();

	for (;A_iterator != A_end; ++A_iterator)
		OS_Array_xor(*A_iterator,right);
	return result;
};

/*! \brief Matrix multiplication A*B
	Perform the algebra operation A*B, where A and B are two matrixes.
	If \f$ A = \{a_{i,j}}_{i,j}\f$ and  \f$ B = \{b_{i,j}}_{i,j}\f$, then it returns the matrix \f$ \{ \sum_k a_{i,k} \cdot b_{k,j}\}_{i,j} \f$.
*/
template<class A,class B>
inline B OS_Array_Matrix_Multiply(const A& left, const B& right) {
	OS_STATIC_CHECK(A::Dim == 2*B::Dim);
	typedef typename A::const_iterator A_const_iterator_type;
	typedef typename B::iterator B_iterator_type;

	B result(right.size_list());
	A_const_iterator_type A_const_iterator		= A.begin();
	A_const_iterator_type A_const_iterator_end	= A.end();
	B_iterator			  result_iterator		= result.begin();

	for (;A_const_iterator != A_const_iterator_end; ++A_const_iterator, ++result_iterator)
		*result_iterator = (*A_const_iterator)* (right);
	return result;
};


/*! \brief partial scalar multiplication A*b
	Perform a partial operation of scalar multiplication.
	If \f$ A = \{a_i}_i\f$ then it corresponds to \f$ \{a_i \cdot b\}_{I \le i < J} \f$.
*/
template<class Iterator, typename value_type>
inline void OS_Array_multiply(Iterator left, Iterator right, value_type& from)
{
	for (;left != right; ++left)
		*left *= from;
};

/*! \brief scalar multiplication A*b
	Perform a scalar multiplication.
	If \f$ A = \{a_i}_i\f$ then it corresponds to \f$ \{a_i \cdot b\}_{i} \f$.
*/
template<class A>
inline A& OS_Array_multiply(A& left, const typename A::base_value_type& from)
{
	typedef typename A::size_type size_type;
	typedef typename A::base_value_type base_value_type;
	base_value_type * begin = left.first_;
	base_value_type * end = left.first_+left.capacity();
	for (;begin != end; ++begin)
		*begin *= from;
	return left;
};



/*! \brief scalar division A/b
	Perform a scalar multiplication.
	If \f$ A = \{a_i}_i\f$ then it corresponds to \f$ \{ \frac{a_i}{b}\}_{i} \f$.
*/
template<class A>
inline A& OS_Array_divide(A& left, const typename A::base_value_type from)
{
	typedef typename A::size_type size_type;
	typedef typename A::base_value_type base_value_type;
	base_value_type * begin = left.first_;
	base_value_type * end = left.first_+left.capacity();
	for (;begin != end; ++begin)
		*begin /= from;
	return left;
};
/*! \brief scalar division A/b
	Perform a scalar multiplication.
	If \f$ A = \{a_i}_i\f$ then it corresponds to \f$ \{ \frac{a_i}{b}\}_{i} \f$.
*/
template<class A>
inline A& OS_Array_divide(const typename A::base_value_type from, A& left)
{
	typedef typename A::size_type size_type;
	typedef typename A::base_value_type base_value_type;
	base_value_type * begin = left.first_;
	base_value_type * end = left.first_+left.capacity();
	A::base_value_type local;
	for (;begin != end; ++begin) {
		local = *begin;
		if (local == 0.)
			*begin = from / THRESHOLD;
		else 
			*begin = from / local;
	};
	return left;
};

template<class A, class B>
inline A& OS_Array_external_multiply(A& left, const B& from)
{
	typedef typename A::value_type value_type;
	typedef typename A::iterator iterator;
	OS_STATIC_CHECK(A::Dim-1 == B::Dim);
	
	iterator beg = left.begin();
	iterator end = left.end();
	for (;beg != end; ++beg) {
		OS_Array_xor( (*beg), from); 
	};
	return left;
};

template<class A, class B>
inline A& OS_Array_external_divide(A& left, const B& from)
{
	return OS_Array_external_multiply(left,1./ from);
};
/*! \brief multiplication by -1.
	If \f$ A = \{a_i}_i\f$ then it returns to \f$ \{ -a_i \}_{i} \f$.
*/
template<class A>
inline A OS_Array_negate(A& left)
{	
	typedef typename A::size_type size_type;
	typedef typename A::base_value_type base_value_type;
	A result( (const A&) left);
	base_value_type * begin = result.first_;
	base_value_type * end = result.first_+result.capacity();
	for (;begin != end; ++begin)
		*begin *= -1.;
	return result;
};
/*! \brief inversion of a vector
	If \f$ A = \{a_i}_i\f$ then it returns to \f$ \{ \frac{1}{a_i}\}_{i} \f$.
*/
template<class A>
inline A& OS_Array_not(A& left)
{	
	typedef typename A::size_type size_type;
	typedef typename A::base_value_type base_value_type;
	base_value_type * begin = left.first_;
	base_value_type * end = left.first_+left.capacity();
	for (;begin != end; ++begin)
		*begin = 1/(*begin);
	return left;
};
/*! \brief power of a vector
	If \f$ A = \{a_i}_i\f$ then it returns to \f$ \{ (a_i)^b \}_{i} \f$.
*/
template<class A, typename D>
inline A& OS_Array_pow(A& left, const D& right)
{	
	typedef typename A::base_value_type base_value_type;
	base_value_type * begin = left.first_;
	base_value_type * end = left.first_+left.capacity();
	for (;begin != end; ++begin)
		*begin = std::pow((base_value_type&) *begin,(D&) right);
	return left;
};
/*! \brief exponential of a vector
	If \f$ A = \{a_i}_i\f$ then it returns to \f$ \{ \exp(a_i) \}_{i} \f$.
*/
template<class A>
inline A& OS_Array_exp(A& left)
{	
	typedef typename A::base_value_type base_value_type;
	base_value_type * begin = left.first_;
	base_value_type * end = left.first_+left.capacity();
	for (;begin != end; ++begin)
		*begin = std::exp(*begin);
	return left;
};
/*! \brief square root of a vector
	If \f$ A = \{a_i}_i\f$ then it returns to \f$ \{ \sqrt{a_i} \}_{i} \f$.
*/
template<class A>
inline A& OS_Array_sqrt(A& left)
{	
	typedef typename A::base_value_type base_value_type;
	base_value_type * begin = left.first_;
	base_value_type * end = left.first_+left.capacity();
	for (;begin != end; ++begin)
		*begin = std::sqrt(*begin);
	return left;
};
/*! \brief logarithm of a vector
	If \f$ A = \{a_i}_i\f$ then it returns to \f$ \{ \ln(a_i) \}_{i} \f$.
*/
template<class A>
inline A& OS_Array_ln(A& left)
{	
	typedef typename A::base_value_type base_value_type;
	base_value_type * begin = left.first_;
	base_value_type * end = left.first_+left.capacity();
	for (;begin != end; ++begin)
		*begin = std::log(*begin);
	return left;
};
/*! \brief maximum of two vectors 
	If \f$ A = \{a_i}_i\f$ and \f$ B = \{b_i}_i\f$ then it returns to \f$ \{ \sup(a_i,b_i) \}_{i} \f$.
*/
template<class A, class B>
inline A& OS_Array_max(A& left, const B& right)
{	
	OS_size cap = left.capacity();
	OS_DYNAMIC_CHECK( cap ==right.capacity(), "OS_Array_Copy : this function suppose that the destination vector is correctly sized" );

	typedef typename A::base_value_type base_value_type;
	typedef typename B::base_value_type base_value_type_B;
	base_value_type * begin = left.first_;
	base_value_type_B * begin_B = right.first_;
	base_value_type * end = left.first_+cap;
	for (;begin != end; ++begin, ++begin_B)
		*begin = std::max(*begin,*begin_B);
	return left;
};

/*! \brief vectors  affectation
	If \f$ A = \{a_i}_i\f$ then it fills it with b: \f$ \{ a_i = b \}_{i} \f$.
*/
template<class A, class fundamental = boost::is_fundamental<A>::type >
struct OS_Array_assign
{
	void operator()(A& left,const typename A::base_value_type from)
	{	
		typedef typename A::size_type size_type;
		typedef typename A::base_value_type base_value_type;
		base_value_type * begin = left.first_;
		base_value_type * end = left.first_+left.capacity();
		for (;begin != end; ++begin)
			*begin = from;
	};
};

template<class A>
struct OS_Array_assign<A,OS_true>
{
	void operator()(A& left, A from)
	{	
		left=from;
	};
};

/*! \brief vector copy
*/
template<class A, class B >
	A& OS_Array_Copy(A& dest,const B & source)
	{	
		OS_size cap = dest.capacity();
		OS_DYNAMIC_CHECK( cap ==source.capacity(), "OS_Array_Copy : this function suppose that the destination vector is correctly sized" );
		typedef A::base_value_type base_value_type;
		memcpy(dest.first_,source.first_,cap*sizeof(base_value_type));
		return dest;
	};

/*! \brief $\ell^2$ norm of vectors
	If \f$ A = \{a_i}_i\f$ then it returns b: \f$ \sum_i a_i^a_i \f$.
*/
template<class A>
inline typename A::value_type OS_norm22(const A& from)
{
	typedef typename A::value_type value_type;
	typedef typename A::const_iterator const_iterator;
	typedef typename value_type::iterator iterator;
	typedef typename A::indice_list_type indice_list_type;

	A temp(from^from);
	indice_list_type local_size(from.size_list());
	local_size.pop_back();
	value_type result(local_size);
	result.assign(0.);
	const_iterator beg = temp.begin();
	const_iterator end = temp.end();
	for (;beg != end; ++beg) {
		result += *beg; 
	};
	return result;
};
/*! \brief square root of $\ell^2$ norm of vectors
	If \f$ A = \{a_i}_i\f$ then it returns b: \f$ \sqrt{ \sum_i a_i^a_i }\f$.
*/
template<class A>
inline typename A::value_type OS_norm2(const A& from)
{
	return sqrt(OS_norm22(from));
};

/*! \brief binary boolean equality operator
	If \f$ A = \{a_i}_i\f$ and \f$ B = \{B_i}_i\f$ then it returns true if \f$ {a_i == b_i}\f$ for all i.
*/
template<class A,class B>
bool OS_Array_equal(const A& left, const B& right)
{
	OS_STATIC_CHECK(A::Dim>0 && B::Dim>0);
	typedef typename A::const_iterator const_iterator_A;
	typedef typename A::const_iterator const_iterator_A;
	typedef typename B::const_iterator const_iterator_B;
	for (int i = 1;i<= A::Dim;++i)
		if(left.size(i) != right.size(i)) 
			return false;
	const_iterator_A it_A = left.begin();
	const_iterator_B it_B = right.begin();
	const_iterator_A it_end = left.end();
	for (;it_A != it_end; ++it_A, ++it_B)
		if ( (*it_A) != (*it_B) ) 
			return false;
	return true;
};
/*! \brief scalar product of vectors
	If \f$ A = \{a_i}_i\f$ and \f$ B = \{B_i}_i\f$ then it returns true if \f$ \sum a_i* b_i \f$ for all i.
*/
template<class A,class B>
inline OS_double OS_Array_scalar_product(const A & left, const B & right) {
		typedef typename A::size_type size_type;
		typedef typename A::base_value_type base_value_type;
		for (int i = 1;i<= A::Dim;++i)
			OS_DYNAMIC_CHECK(left.size(i) == right.size(i),"scalar product of inconsistent OS_array"); // not enough but improbable
		OS_double product = 0;
		base_value_type * from1_begin = left.first_;
		base_value_type * from2_begin = right.first_;
		base_value_type * from1_end = from1_begin + left.capacity();
		for (;from1_begin != from1_end; ++from1_begin, ++from2_begin)
			product += (*from1_begin) * (* from2_begin);
		return product;
};

/*! \brief binary boolean \f$ \le \f$ operator
	If \f$ A = \{a_i}_i\f$ and \f$ B = \{B_i}_i\f$ then it returns true if \f$ {a_i \le b_i}\f$ for all i.
*/
template<class A,class B>
struct OS_Array_less 
{
	OS_Array_less (const A& left, const B& right) : result_(true)
	{
		OS_STATIC_CHECK(A::Dim>0 && B::Dim>0);
		typedef typename A::const_iterator const_iterator_A;
		typedef typename B::const_iterator const_iterator_B;
		OS_DYNAMIC_CHECK(left.size_list() == right.size_list(),"comparing inconsistent OS_array");
		const_iterator_A it_A = left.begin();
		const_iterator_B it_B = right.begin();
		const_iterator_A it_end = left.end();
		for (;it_A != it_end; ++it_A, ++it_B)
			if ( !( (*it_A) < (*it_B) ) ) {
				result_ = false;
				break;
			};
	};
	inline operator bool() {return result_;};
	bool result_;
};
/*! \brief unary boolean operator \f$ \le \f$ operator for vectors
	If \f$ A = \{a_i}_i\f$ then it returns true if \f$ {a_i \le b}\f$ for all i.
*/
template<class A>
struct OS_Array_less<A, typename A::base_value_type> {
	inline operator bool() {return result_;};
	OS_Array_less (const A& left, typename A::base_value_type right) : result_(true)
	{
		OS_STATIC_CHECK(A::Dim>0);
		typedef typename A::const_iterator const_iterator_A;
		const_iterator_A it_A = left.begin();
		const_iterator_A it_end = left.end();
		for (;it_A != it_end; ++it_A)
			if ( !( (*it_A) < right ) ) {
				result_ = false;
				break;
			};
	};
	bool result_;
};

/*! \brief unary boolean operator \f$ \le \f$ operator for vectors
	If \f$ A = \{a_i}_i\f$ then it returns true if \f$ {a_i \ge b}\f$ for all i.
*/
template<class A>
struct OS_Array_less<typename A::base_value_type, A> {
	inline operator bool() {return result_;};
	OS_Array_less (typename A::base_value_type left, const A& right) : result_(true)
	{
		OS_STATIC_CHECK(A::Dim>0);
		typedef typename A::const_iterator const_iterator_A;
		const_iterator_A it_A = right.begin();
		const_iterator_A it_end = right.end();
		for (;it_A != it_end; ++it_A)
			if ( !( (*it_A) > left ) ) {
				result_ = false;
				break;
			};
	};
	bool result_;
};

////////////////////////////////////////////


template< template< typename> class T, OS_size Dim, class base_value_type>
struct OS_Multi_Array {
	typedef T<typename OS_Multi_Array<T,Dim-1,base_value_type>::type> type;
};
template<template< typename> class T,class base_value_type>
struct OS_Multi_Array<T,1,base_value_type> {
	typedef T<base_value_type> type;
};
template< template< typename> class T,class base_value_type>
struct OS_Multi_Array<T,0,base_value_type> {
	typedef base_value_type type;
};


template<typename T> 
struct OS_data_array;

template <class T>
struct OS_array_type_traits
{
		typedef T														base_value_type;
		typedef T*														ptr_type;
		typedef OS_int													indice_type;
		typedef OS_size													size_type;
		typedef std::vector<size_type>									indice_list_type;
		typedef std::vector<size_type>									capacity_list_type;
		typedef OS_data_array<T>										data_type;
		typedef smart_ptr<data_type>									smart_ptr_type;
};
/*! \brief Basic class of contiguous set of datas for Multidimensional array
*/
template<typename T> 
struct OS_data_array : public OS_array_type_traits<T>
{

		OS_data_array(const OS_data_array &from, const size_type dim , const ptr_type first ) : size_(dim), capacity_(dim), dim_(dim)
		{
			std::copy(from.size_.begin(),from.size_.begin()+dim_,size_.begin());
			std::copy(from.capacity_.begin(),from.capacity_.begin()+dim_,capacity_.begin());
			allocate();
			OS_DYNAMIC_CHECK( (first-from.first_) + capacity() <= from.capacity(), "OS_data_array : bad assignement parameters");
			if (first_)
				memcpy(first_, first, capacity()*sizeof(base_value_type));
		};
		OS_data_array(const OS_data_array &from) : dim_(from.dim_), size_(from.size_), capacity_(from.capacity_)
		{
			allocate();
			if (first_)
				memcpy(first_, from.first_, capacity()*sizeof(base_value_type));
		};

		OS_data_array(): first_(0), dim_(0) {
		};
		~OS_data_array(){if (first_) delete first_;};
		explicit OS_data_array(size_type nb) : dim_(1), size_(1), capacity_(1)
		{
			size_[0] = nb;
			init_capacity();
			allocate();
		};
		explicit OS_data_array(size_type dim1, size_type dim2) :dim_(2), size_(2), capacity_(2)
		{
			size_[0] = dim2;
			size_[1] = dim1;
			init_capacity();
			allocate();
		};
		explicit OS_data_array(size_type dim1, size_type dim2, size_type dim3) :dim_(3), size_(3), capacity_(3) {
			size_[0] = dim3;
			size_[1] = dim2;
			size_[2] = dim1;
			init_capacity();
			allocate();
		};

		OS_data_array( const indice_list_type& from) : size_(from), dim_(from.size()), capacity_(from.size()) { 
			init_capacity();
			allocate();
		};
		OS_data_array& operator =(const OS_data_array &from) {
			if (size_ !=	from.size_list()) {
				dim_ = from.dim_;
				size_		= from.size_list();
				capacity_	= from.capacity_list();
				deallocate();
				allocate();
			};
			if (first_)
				memcpy(first_, from.first_, capacity()*sizeof(base_value_type));
			return *this;
		};

		void resize(const indice_type & dim1) {
			if (dim_ != 1) {
				dim_ = 1;
				size_.resize(1);
				capacity_.resize(1);
			};
			if (size_[0] != dim1) {
				size_[0] = dim1;
				reallocate();
			};
		};
		void resize(const indice_type & dim1,const indice_type & dim2) {
			if (dim_ != 2) {
				dim_ = 2;
				size_.resize(dim_);
				capacity_.resize(dim_);
			};
			if (size_[0] != dim2 | size_[1] != dim1) {
				size_[0] = dim2;size_[1] = dim1;
				reallocate();
			};
		};
		void resize(const indice_type & dim1,const indice_type & dim2,const indice_type & dim3) {
			if (dim_ != 3) {
				dim_ = 3;
				size_.resize(dim_);
				capacity_.resize(dim_);
			};
			if (size_[0] != dim3 | size_[1] != dim2 | size_[2] != dim1) {
				size_[0] = dim3;size_[1] = dim2;size_[2] = dim1;
				reallocate();
			};
		};

		void resize(const indice_list_type & indx) {
			if (dim_ != indx.size()) {
				dim_ = indx.size();
				size_.resize(dim_);
				capacity_.resize(dim_);
			};
			if (size_ != indx) {
				dim_ = indx.size();
				size_ = indx;
				capacity_.resize(dim_);
				reallocate();
			};
		};
		inline size_type dim(){return dim_;};
		inline const indice_list_type & size_list() const				{return size_;};
		inline const capacity_list_type & capacity_list() const			{return capacity_;};
		inline size_type size() const									{return size_[dim_-1];} ;
		inline size_type size(size_type which_dim) const {
			if (dim_ == 0)
				return 0;
			if (which_dim == 0)
				return 1;
			return size_[which_dim-1];
		};
		inline size_type capacity(size_type which_dim) const {
			if (dim_ == 0)
				return 0;
			if (which_dim == 0)
				return 1;
			return capacity_[which_dim-1];
		};
		inline size_type capacity() const									{return capacity_[dim_-1];};
		void reallocate() {
			deallocate();
			init_capacity();
			allocate();
		};

private :

		inline void init_capacity() {
		  capacity_[0] = size_[0];
		  for (OS_int i = 1 ;i < dim_ ; ++i)
			 capacity_[i] = capacity_[i-1]*size_[i];
		};
		inline void deallocate() {
			if (first_) delete first_;
		};
		inline void allocate() {
			size_type cap(0);
			if (dim_) {
				cap = capacity();
			};
			if (cap) 
				first_ = new T[cap]; 
			else 
				first_ = 0;
		};
public :
		ptr_type				first_;
		indice_list_type		size_;
		capacity_list_type		capacity_;
		size_type dim_;
};

template<typename A, std::size_t B >
class OS_array;



/*! \brief Basic class for references of Multidimensional array
*/
template<typename T,OS_size Dim = T::Dim>
class OS_data_array_ref  :
	public OS_array_type_traits<typename T::base_value_type>
{
public :
	enum :OS_size{Dim = Dim};
	OS_STATIC_CHECK(Dim>0);
	typedef typename T												cast_type;
	typedef typename cast_type::value_type							value_type;
	typedef typename OS_data_array_ref								type;
//////////////////////////////// reference
	template <size_type dim>
	struct reference_helper {
		OS_STATIC_CHECK(dim>0);
		typedef typename OS_data_array_ref<value_type> reference;
		typedef typename const OS_data_array_ref<value_type>& const_reference;
	};
	template <>
	struct reference_helper<1> {
		typedef typename base_value_type& reference;
		typedef typename const base_value_type& const_reference;
	};
	typedef  typename cast_type::reference					reference;
	typedef  typename cast_type::const_reference			const_reference;
///////////////////////////

	typedef typename cast_type::indice_list_type							indice_list_type;
	typedef typename cast_type::capacity_list_type							capacity_list_type;
	typedef typename cast_type::iterator									iterator;
	typedef typename cast_type::const_iterator								const_iterator;
	
	OS_data_array_ref() : ref_(0),first_(0){};
	explicit OS_data_array_ref(cast_type &  from) : ref_(from.data_),first_(from.first_){};

	explicit OS_data_array_ref(smart_ptr_type ref, ptr_type value) : ref_(ref),first_(value){};
	type & operator = (const cast_type & from){
		OS_DYNAMIC_CHECK(capacity() == from.capacity(), "OS_data_array_ref : bad assignement parameters");
		memcpy(first_, from.first_, capacity()*sizeof(base_value_type));
		return *this;
	};
	type & operator = (const type & from){
		OS_DYNAMIC_CHECK(capacity() == from.capacity(), "OS_data_array_ref : bad assignement parameters");
		memcpy(first_, from.first_, capacity()*sizeof(base_value_type));
		return *this;
	};
	inline const indice_list_type & size_list()	const			{if (ref_) return ref_->size_list(); else return indice_list_type();};
	inline const capacity_list_type & capacity_list() const		{if (ref_) return ref_->capacity_list();else return capacity_list_type();};
	inline size_type size() const								{if (ref_) return ref_->size(Dim);else return 0;};
	inline size_type size(size_type which_dim)	const			{if (ref_) return ref_->size(which_dim);else return 0;};
	inline size_type capacity(size_type which_dim)	const		{if (ref_) return ref_->capacity(which_dim);else return 0;};
	inline size_type capacity() const							{if (ref_) return ref_->capacity(Dim);else return 0;};

	inline type& operator += ( const type &  from){
		return OS_Array_plus(*this,from);
	};
	inline type& operator += ( const cast_type &  from){
		return OS_Array_plus(*this,from);
	};
	inline type& operator -= ( const type &  from){
		return OS_Array_minus(*this,from);
	};
	inline type& operator -= ( const cast_type &  from) {
		return OS_Array_minus(*this,from);
	};
    inline type& operator*=(const base_value_type & from)  {
		return OS_Array_multiply(*this,from);
     };
    type& operator/=(const base_value_type & from)  {
		return OS_Array_divide(*this,from);
     };
	inline type& operator ^= ( const type &  from){
		return 	OS_Array_xor(*this,from);
	};
	inline type& operator ^= ( const cast_type &  from){
		return 	OS_Array_xor(*this,from);
	};
    inline cast_type operator -()  const {
		return -cast_type(*this);
     };
   inline cast_type operator !()  const {
		return -cast_type(*this);
     };
	inline void assign (const base_value_type & from) {
		OS_Array_assign<type>()(*this,from);
	}
	inline bool operator < ( const type &  from) const {
		return OS_Array_less<type,type>(*this,from);
	};
	inline bool operator == ( const type &  from) const {
		return OS_Array_equal(*this,from);
	};
	inline bool operator == ( const cast_type &  from) const {
		return OS_Array_equal(*this,from);
	};
	inline bool operator !=(const type & from) const {
		return ! (*this == from);
	};
	inline type& operator <<= (std::size_t how) {
		cast_type(ref_,first_)<<= how;
		return *this;
	};
	inline type& operator >>= (std::size_t how) {
		cast_type(ref_,first_)>>= how;
		return *this;
	};

	reference operator[](indice_type n) {
		return cast_type(ref_,first_)[n];
	};
	const_reference operator[](indice_type n) const {
		return cast_type(ref_,first_)[n];
	};

	inline typename cast_type::iterator begin() {
		return cast_type(ref_,first_).begin();
	};
	inline typename cast_type::iterator end() {
		return cast_type(ref_,first_).end();
	};
	inline typename cast_type::const_iterator begin() const {
		return cast_type(ref_,first_).begin();
	};
	inline typename cast_type::const_iterator end() const {
		return cast_type(ref_,first_).end();
	};


	smart_ptr_type ref_;
	mutable ptr_type	first_;
};

/*! \brief right addition operator of a reference of Multidimensional array and a Multidimensional Array
*/
template<typename T,OS_size D>
		inline typename OS_data_array_ref<T,D>::cast_type operator +(const OS_data_array_ref<T,D> & from1, const typename OS_data_array_ref<T,D>::cast_type & from2) {
		typedef OS_data_array_ref<T,D>::cast_type A;
		A temp(from1);
		::OS_Array_plus(temp,from2);
		return temp;
};
/*! \brief left addition operator of a Multidimensional Array and a Multidimensional array
*/
template<typename T,OS_size D>
		inline typename OS_data_array_ref<T,D>::cast_type operator +(const typename OS_data_array_ref<T,D>::cast_type & from2,const OS_data_array_ref<T,D> & from1) {
		typedef OS_data_array_ref<T,D>::cast_type A;
		A temp(from1);
		::OS_Array_plus(temp,from2);
		return temp;
};
/*! \brief addition operator of a Multidimensional Array references.
*/
template<typename T,OS_size D>
		inline typename OS_data_array_ref<T,D>::cast_type operator +(const OS_data_array_ref<T,D> & from1, const OS_data_array_ref<T,D> & from2) {
		typedef OS_data_array_ref<T,D>::cast_type A;
		A temp(from1);
		::OS_Array_plus(temp,from2);
		return temp;
};
/*! \brief minus operator of Multidimensional Array references.
*/
template<typename T,OS_size D>
		inline typename OS_data_array_ref<T,D>::cast_type operator -(const OS_data_array_ref<T,D> & from1, const OS_data_array_ref<T,D> & from2) {
		typedef OS_data_array_ref<T,D>::cast_type A;
		A temp(from1);
		OS_Array_minus(temp,from2);
		return temp;
};
/*! \brief right minus operator of a Multidimensional Array reference and a Multidimensional Array.
*/
template<typename T,OS_size D>
		inline typename OS_data_array_ref<T,D>::cast_type operator -(const OS_data_array_ref<T,D> & from1, const typename OS_data_array_ref<T,D>::cast_type & from2) {
		typedef OS_data_array_ref<T,D>::cast_type A;
		A temp(from1);
		::OS_Array_minus(temp,from2);
		return temp;
};
/*! \brief left minus operator of a Multidimensional Array and a Multidimensional Array reference.
*/
template<typename T,OS_size D>
		inline typename OS_data_array_ref<T,D>::cast_type operator -(const typename OS_data_array_ref<T,D>::cast_type & from2,const OS_data_array_ref<T,D> & from1) {
		typedef OS_data_array_ref<T,D>::cast_type A;
		A temp(from1);
		::OS_Array_minus(temp,from2);
		return temp;
};

/*! \brief external product operator of Multidimensional Array reference.
*/
template<typename T, std::size_t D>
inline typename OS_data_array_ref<T,D>::cast_type operator ^(const OS_data_array_ref<T,D> & left, const OS_data_array_ref<T,D> & right) {
		OS_data_array_ref<T,D>::cast_type temp(left);
		OS_Array_xor(temp,right);
		return temp;
};
template<typename T, std::size_t D>
inline typename OS_data_array_ref<T,D>::cast_type operator ^(const OS_data_array_ref<T,D> & left, const typename OS_data_array_ref<T,D>::cast_type & right) {
		OS_data_array_ref<T,D>::cast_type temp(left);
		OS_Array_xor(temp,right);
		return temp;
};
/*! \brief scalar product operator of Multidimensional Array reference.
*/
template<typename T, std::size_t D>
inline OS_double operator *(const OS_data_array_ref<T,D> & left, const OS_data_array_ref<T,D> & right) {
		return OS_Array_scalar_product(left,right);
};
/*! \brief matrix multiplication of Multidimensional Array reference.
*/
template<typename T, std::size_t D>
inline typename OS_data_array_ref<T,D>::cast_type operator *(const OS_data_array_ref<T,2*D> & left, const OS_data_array_ref<T,D> & right) {
		return OS_Array_Matrix_Multiply(left,right);
};
/*! \brief division of Multidimensional Array reference.
*/
template<typename T,OS_size D>
inline typename OS_data_array_ref<T,D>::cast_type operator /(const OS_data_array_ref<T,D> & from1, const typename OS_data_array_ref<T,D>::cast_type & from2) {
		typedef OS_data_array_ref<T,D>::cast_type A;
		A temp(from1);
		OS_Array_xor_not(temp,from2);
		return temp;
};
template<typename T,OS_size D>
inline typename OS_data_array_ref<T,D>::cast_type operator /(const typename OS_data_array_ref<T,D>::cast_type & from1, const OS_data_array_ref<T,D> & from2) {
		typedef OS_data_array_ref<T,D>::cast_type A;
		A temp(from1);
		OS_Array_xor_not(temp,from2);
		return temp;
};
template<typename T,OS_size D>
inline typename OS_data_array_ref<T,D>::cast_type operator /(const OS_data_array_ref<T,D> & from1, const typename OS_data_array_ref<T,D>::base_value_type & from2) {
		typedef OS_data_array_ref<T,D>::cast_type A;
		A temp(from1);
		OS_Array_divide(temp,from2);
		return temp;
};
template<typename T,OS_size D>
inline typename OS_data_array_ref<T,D>::cast_type operator /(const typename OS_data_array_ref<T,D>::base_value_type & from2, const OS_data_array_ref<T,D> & from1) {
		typedef OS_data_array_ref<T,D>::cast_type A;
		A temp(from1);
		OS_Array_divide(temp,from2);
		return temp;
};
template<typename T,OS_size D>
inline typename OS_data_array_ref<T,D>::cast_type operator *(const OS_data_array_ref<T,D> & from1, const typename OS_data_array_ref<T,D>::base_value_type & from2) {
		typedef OS_data_array_ref<T,D>::cast_type A;
		A temp(from1);
		OS_Array_multiply(temp,from2);
		return temp;
};
template<typename T,OS_size D>
inline typename OS_data_array_ref<T,D>::cast_type operator *(const typename OS_data_array_ref<T,D>::base_value_type & from2, const OS_data_array_ref<T,D> & from1) {
		typedef OS_data_array_ref<T,D>::cast_type A;
		A temp(from1);
		OS_Array_multiply(temp,from2);
		return temp;
};
template<typename T,OS_size D>
inline typename OS_data_array_ref<T,D>::cast_type operator <<(const OS_data_array_ref<T,D> & from, std::size_t n) {
		typedef OS_data_array_ref<T,D>::cast_type A;
		A temp(from);
		return temp<<=n;
};

template<typename T,OS_size D>
inline typename OS_data_array_ref<T,D>::cast_type operator >>(const OS_data_array_ref<T,D> & from, std::size_t n) {
		typedef OS_data_array_ref<T,D>::cast_type A;
		A temp(from);
		return temp>>=n;
};





/*! \brief serialization operator for Multidimensional Array reference.
*/
template<typename T,OS_size D>
std::ostream &operator<<(std::ostream &out, const OS_data_array_ref<T,D> &v) {
	typedef typename OS_data_array_ref<T,D>::cast_type cast_type;
	cast_type breakpoint(v.ref_,v.first_);
	return out << breakpoint;
};


/*! \brief Iterator for Multidimensional Array.
*/
template<typename T,OS_size Dim = T::Dim>
class OS_array_iterator 
	:public OS_array_type_traits<typename T::base_value_type>
	, public boost::random_access_iterator_helper<
     OS_array_iterator<T, Dim>
	 , typename T
	 , std::ptrdiff_t
	 , typename T::ptr_type
	 , typename OS_data_array_ref<T,Dim>
	>
{
public :
	enum :OS_size{Dim = Dim};
	OS_STATIC_CHECK(Dim>=0);
	typedef OS_array_iterator										type;
	typedef typename T												value_type;
	typedef typename OS_data_array_ref<T,Dim>						reference;

	OS_array_iterator (const OS_array_iterator & from) : ref_(from.ref_), it_value_(from.it_value_), capacity_(from.capacity_)
	{}; 
	OS_array_iterator() : it_value_(0), capacity_(0) {};
	explicit OS_array_iterator(smart_ptr_type ref, ptr_type value) : ref_(ref),it_value_(value), capacity_(ref_->capacity(Dim)){};

	inline reference operator*(){
		return reference(ref_,it_value_);
	};
	inline OS_array_iterator& operator +=(size_type n) {
		it_value_ += n*capacity_;
		return *this;
	};
	inline OS_array_iterator& operator -=(size_type n) {
		it_value_ -= n*capacity_;
		return *this;
	};
	inline OS_array_iterator operator -(size_type n) const {
		OS_array_iterator result(*this);
		result.it_value_ -= n*capacity_;
		return result;
	};
	inline OS_array_iterator operator +(size_type n) const {
		OS_array_iterator result(*this);
		result.it_value_ += n*capacity_;
		return result;
	};
	inline OS_array_iterator& operator ++() {
		it_value_ += capacity_;
		return *this;
	};
	inline OS_array_iterator& operator --() {
		it_value_ -= capacity_;
		return *this;
	};
	inline bool operator ==(const OS_array_iterator & from) const {
		return (it_value_ == from.it_value_ && ref_ == from.ref_);
	};
	inline bool operator<(const OS_array_iterator & x) const {
		return it_value_ < x.it_value_;
	};
	smart_ptr_type		ref_;
	size_type capacity_;
	mutable ptr_type	it_value_;
};

template<typename A , std::size_t D>
std::ptrdiff_t operator-(const OS_array_iterator<A,D>& x, 
						 const OS_array_iterator<A,D>& y) {
							 OS_DYNAMIC_CHECK(x.ref_ == y.ref_,"std::ptrdiff_t operator-(const OS_array_iterator<A,D> : Bad paramaters ");  // I have to think a little bit
    return (std::ptrdiff_t)(x.it_value_ - y.it_value_) / x.ref_->capacity(D) ; 
 };

/*! \brief Multidimensional Array.
*/
template<typename A, std::size_t Dim = OS_type_helper<A>::Dim +1 >
class OS_array : public OS_array_type_traits<typename OS_type_helper<A>::type::base_value_type>
{
public :

//	enum :OS_size{Dim = D};
	enum :OS_size{Dim = OS_type_helper<A>::type::Dim +1 };
	typedef OS_array													type;
	typedef OS_array_type_traits										Base;
	static std::string get_name() {return " OS_Array ";};
///////////////////////////value_type

	template <size_type D>
	struct value_type_helper {
		typedef A type;
	};
	template <>
	struct value_type_helper<1> {
		typedef base_value_type type;
	};
	typedef typename value_type_helper<Dim>::type					value_type;
/////////////// iterator
	template <size_type dim>
	struct iterator_helper {
		OS_STATIC_CHECK(dim>0);
		typedef typename OS_array_iterator<value_type> type;
	};
	template <>
	struct iterator_helper<1> {
		typedef base_value_type* type;
	};
	typedef  typename iterator_helper<Dim>::type					iterator;
	typedef  typename iterator_helper<Dim>::type					const_iterator;

//////////////////////////////// reference
	template <size_type dim>
	struct reference_helper {
		OS_STATIC_CHECK(dim>0);
		typedef typename OS_data_array_ref<value_type>			type;
		typedef const typename OS_data_array_ref<value_type>	const_type;
	};
	template <>
	struct reference_helper<1> {
		typedef const base_value_type& const_type;
		typedef base_value_type& type;
	};

	typedef  typename reference_helper<Dim>::type					reference;
	typedef  typename reference_helper<Dim>::const_type				const_reference;


	explicit OS_array(size_type nb) : data_(new data_type(nb)) {
		OS_STATIC_CHECK(Dim==1);
		first_ = data_->first_;
	};
	explicit OS_array(size_type dim1, size_type dim2) : data_(new data_type(dim1,dim2)) {
		OS_STATIC_CHECK(Dim==2);
		first_ = data_->first_;
	};
	explicit OS_array(size_type dim1, size_type dim2, size_type dim3) : data_(new data_type(dim1,dim2,dim3)) {
		OS_STATIC_CHECK(Dim==3);
		first_ = data_->first_;
	};
	explicit OS_array(const indice_list_type& indx) : data_(new data_type(indx)) {first_ = data_->first_;};
	inline void resize(size_type nb) {data_->resize(nb); first_ = data_->first_;};
	inline void resize(size_type dim1, size_type dim2) {data_->resize(dim1,dim2); first_ = data_->first_;};
	inline void resize(size_type dim1, size_type dim2, size_type dim3) {data_->resize(dim1,dim2,dim3); first_ = data_->first_;};
	inline void resize(const indice_list_type& indx) {data_->resize(indx); first_ = data_->first_;};

	explicit OS_array(smart_ptr_type data, ptr_type first) : data_(data), first_(first){};

	OS_array() : data_(new data_type){first_ = data_->first_;};
	virtual ~OS_array(){};
	OS_array(const OS_array & from) : data_(new data_type(*from.data_)){first_ = data_->first_;};
	OS_array& operator = (const OS_array & from) 
	{
		data_->resize(from.data_->size_list());
		*data_ = *from.data_;
		first_ = data_->first_;
		return *this;
	};

	template <class T>
	OS_array(const OS_data_array_ref<T,Dim> & from) :data_(new data_type(*from.ref_, Dim, from.first_))
	{first_ = data_->first_;};

	inline const indice_list_type & size_list()	const			{if (data_) return data_->size_list(); };
	inline const capacity_list_type & capacity_list() const		{if (data_) return data_->capacity_list();else return capacity_list_type();};
	inline size_type size() const								{if (data_) return data_->size(Dim);else return 0;};
	inline size_type size(size_type which_dim)	const			{if (data_) return data_->size(which_dim);else return 0;};
	inline size_type capacity(size_type which_dim)	const		{if (data_) return data_->capacity(which_dim);else return 0;};
	inline size_type capacity() const							{if (data_) return data_->capacity(Dim);else return 0;};

	inline type& operator += ( const type &  from){
		return OS_Array_plus<type,type>(*this,from); 
	};
	inline type& operator -= ( const type &  from){
		return OS_Array_minus(*this,from); 
	};
    inline type& operator*=(const base_value_type & from)  {
		return OS_Array_multiply(*this,from); 
     };
    inline type& operator/=(const base_value_type & from)  {
		return OS_Array_divide(*this,from); 
     };
    inline type operator -()  const {
		return OS_Array_negate(*this); 
     };
	inline type& operator ^= ( const type &  from){
		return 	OS_Array_xor(*this,from);
	};
    inline type operator !()  const {
		return OS_Array_not(*this); 
     };
	inline void assign (const base_value_type & from) {
		OS_Array_assign<type>()(*this,from);
	}
	inline bool operator == ( const type &  from) const {
		return OS_Array_equal(*this,from);
	};
	inline bool operator !=(const type & from) const {
		return ! (*this == from);
	};

	template<std::size_t Dim>
	iterator begin_helper(){
		return iterator(data_,first_);
	};
	template<>
	iterator begin_helper<1>(){
		return first_;
	};
	inline iterator begin() {
		return begin_helper<Dim>();
	};
	inline iterator end() {
		return begin() + size();
	};
	virtual reference operator[](indice_type n) {
		return *(begin()+n);
	};

	template<std::size_t Dim>
	const_iterator const_begin_helper() const{
		return const_iterator(data_,first_);
	};
	template<>
	const_iterator const_begin_helper<1>() const {
		return first_;
	};
	inline const_iterator begin() const {
		return const_begin_helper<Dim>();
	};
	inline const_iterator end() const {
		return begin() + size();
	};


	smart_ptr_type data_;
	ptr_type first_;
};
template<typename A, std::size_t D>
		inline OS_array<A,D> operator +(const OS_array<A,D> & left, const OS_array<A,D> & right) {
		typedef OS_array<A,D> A;
		A temp(left);
		::OS_Array_plus(temp,right);
		return temp;
};
template<typename A, std::size_t D>
		inline OS_array<A,D> operator -(const OS_array<A,D> & left, const OS_array<A,D> & right) {
		typedef OS_array<A,D> A;
		A temp(left);
		::OS_Array_minus(temp,right);
		return temp;
};

template<typename T,OS_size D>
inline OS_array<T,D> operator /(const OS_array<T,D> & from1, const OS_array<T,D> & from2) {
		typedef OS_array<T,D> A;
		A temp(from1);
		OS_Array_xor_not(temp,from2);
		return temp;
};


template<typename T,OS_size D>
inline OS_array<T,D> operator /(const OS_array<T,D> & from1, const typename OS_array<T,D>::base_value_type & from2) {
		typedef OS_array<T,D> A;
		A temp(from1);
		OS_Array_divide(temp,from2);
		return temp;
};
template<typename T,OS_size D>
inline OS_array<T,D> operator /(const typename OS_array<T,D>::base_value_type & from2, const OS_array<T,D> & from1) {
		typedef OS_array<T,D> A;
		A temp(from1);
		OS_Array_divide(from2,temp);
		return temp;
};
template<typename T,OS_size D>
inline OS_array<T,D> operator *(const typename OS_array<T,D>::base_value_type & from2, const OS_array<T,D> & from1) {
		typedef OS_array<T,D> A;
		A temp(from1);
		OS_Array_multiply(temp,from2);
		return temp;
};
template<typename T,OS_size D>
inline OS_array<T,D> operator *(const OS_array<T,D> & from1, const typename OS_array<T,D>::base_value_type & from2) {
		typedef OS_array<T,D> A;
		A temp(from1);
		OS_Array_multiply(temp,from2);
		return temp;
};
//right tensor contraction
template<typename T, std::size_t D>
inline OS_array<T,D> operator *(const OS_array<T,2*D> & left, const OS_array<T,D> & right) {
		return OS_Array_Matrix_Multiply(left,right);
};

template<typename T,OS_size D, typename A>
inline OS_array<T,D> pow(const OS_array<T,D> & from1, const A& from2) {
		typedef OS_array<T,D> A;
		A temp(from1);
		OS_Array_pow(temp,from2);
		return temp;
};
template<typename T,OS_size D>
inline OS_array<T,D> sqrt(const OS_array<T,D> & from1) {
		typedef OS_array<T,D> A;
		A temp(from1);
		OS_Array_sqrt(temp);
		return temp;
};
template<typename T,OS_size D>
inline OS_array<T,D> exp(const OS_array<T,D> & from1) {
		typedef OS_array<T,D> A;
		A temp(from1);
		OS_Array_exp(temp);
		return temp;
};
template<typename T,OS_size D>
inline OS_array<T,D> ln(const OS_array<T,D> & from1) {
		typedef OS_array<T,D> A;
		A temp(from1);
		OS_Array_ln(temp);
		return temp;
};


template<typename A, std::size_t D>
inline OS_array<A,D> operator ^(const OS_array<A,D> & left, const OS_array<A,D> & right) {
		typedef OS_array<A,D> A;
		A temp(left);
		OS_Array_xor(temp,right);
		return temp;
};


template<typename A, std::size_t D>
inline OS_double operator *(const OS_array<A,D> & left, const OS_array<A,D> & right) {
	return OS_Array_scalar_product(left,right);
};
template<typename A, std::size_t D>
inline bool operator <(const OS_array<A,D> & left, const OS_array<A,D> & right) {
	return ::OS_Array_less<OS_array<A,D>,OS_array<A,D>>(left,right);
};

template<typename A, std::size_t D>
/*! \brief serialization of Multidimensional Array.
*/
std::ostream &operator<<(std::ostream &out, const OS_array<A,D> &v) {
	typedef typename OS_array<A,D>::const_iterator const_iterator;
	const_iterator  begin = v.begin();
	const_iterator	end = v.end();
//	out << "[";
	out << std::endl;
	for (;begin != end; ++begin) {
		out << (*begin);
		if ( begin+1 != end) out << ",";
//		if ( begin+1 != end) {out << "	";} else {out << std::endl;};
	};
//	out << std::endl;
//	out << "]";
  return out;
};


/////////////////////////////////////////////////////:

/*! \brief Multidimensional Vectors.

	It is simply a Multidimensional Array for which a bound check for the operator [] is performed at run time.
*/
template<typename T>
class OS_Vector : 
	public OS_array<T>
{
public :
	typedef	OS_Vector			type;
	typedef OS_array			Base;
	typedef Protected			Boundary_type;

	explicit OS_Vector(size_type nb) : Base (nb) {};
	explicit OS_Vector(size_type dim1, size_type dim2) : Base(dim1,dim2) {};
	explicit OS_Vector(size_type dim1, size_type dim2, size_type dim3) : Base(dim1,dim2,dim3)  {};
	explicit OS_Vector(const indice_list_type& indx) : Base(indx){};
	explicit OS_Vector(smart_ptr_type data, ptr_type first) : Base(data,first){};

	OS_Vector(const Base & from) : Base(from) {};
	OS_Vector(){};
	OS_Vector(const OS_Vector<T>& from)	: Base(from) {};
	OS_Vector& operator = (const type& from)	{
		Base::operator=(from);
		return *this;
	};
	template <class T>
	OS_Vector(const OS_data_array_ref<T,Dim> & from) : Base(from) {};


	virtual typename reference operator[](indice_type n) {
		OS_DYNAMIC_CHECK(n>=0 && n<size(), "indice out of range in array accessor");
		return *(begin()+n);
	};
	virtual typename reference operator[](indice_type n) const {
		OS_DYNAMIC_CHECK(n>=0 && n<size(), "indice out of range in array accessor");
		return *(begin()+n);
	};
};



///////////////////////////

/////////////////////////////

#endif